#!/usr/bin/perl
# 2010/10/05 BedFileで指定した部位のEnrichmentを計算する

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 4 or $ARGV[0] eq '--help'){
	die "Usage : $0 -w [wig file] -b [bed file]\n";
}

my %opt;
getopts("b:w:", \%opt);

my $FILE_bed = $opt{b};
my $FILE_wig = $opt{w};


#---------------------------------------
# method
#---------------------------------------
# 1. wig file を読み込む
# 2. bed file を読み込む
# 3. bed file で指定した部位のenrichを計算する


#---------------------------------------
# read wig file
#---------------------------------------
my %WIG;
#track type=wiggle_0 name=pKu70
#variableStep chrom=chr1 span=30
#50      33.541859
#80      29.758246
#110     26.788264
#140     24.255730
#170     21.654645
#200     19.378484
{
	my $fh;
	if($FILE_wig =~ /\.gz/){
		$fh = IO::File->new("gzip -dc $FILE_wig |") or die "cannot open $FILE_wig: $!";
	}else{
		$fh = IO::File->new($FILE_wig) or die "cannot open $FILE_wig: $!";
	}
	$fh->getline();
	my $chr = '';
	while($_ = $fh->getline()){
		s/\r?\n//;
		if(m/chrom=(\w+) /){
			$chr = $1;
		}else{
			my ($pos, $value) = split /\t/;
			if($value < 0){
				$value = 0;
			}
			$WIG{$chr}{$pos} = $value;
		}
	}
	$fh->close();
}


#---------------------------------------
# read bed file
#---------------------------------------
#chr1    1568383 1568732 LTR.31  0       -       1568383 1568732 0       1       349,    0,
#chr1    2173898 2174235 LTR.35  0       -       2173898 2174235 0       1       337,    0,
#chr1    2474234 2474583 LTR.39  0       -       2474234 2474583 0       1       349,    0,
#chr1    2693885 2694225 LTR.41  0       -       2693885 2694225 0       1       340,    0,
#chr1    2843758 2844091 LTR.42  0       -       2843758 2844091 0       1       333,    0,
my %BED;
my %mask;
{
	my $fh = IO::File->new($FILE_bed) or die "cannot open $FILE_bed: $!";
	while($_ = $fh->getline()){
		s/\r?\n//;
		my ($chr, $start, $end, $t1, $t2, $direction) = split /\t/;
		$BED{$chr}{$start} = [$end, $direction];
		my $mask_start = &convert($start) - 3000;
		my $mask_end = &convert($end) + 3000;
		for(my $i = $mask_start; $i <= $mask_end; $i += 30){
			$mask{$chr}{$i}++;
		}
	}
	$fh->close();
}

### ２つのBED　で示した位置が近すぎる場合は除く (<3kb以内)
my %BED_far;
{
	foreach my $chr(keys %BED){
		foreach my $start(keys %{$BED{$chr}}){
			my ($end, $direction) = @{$BED{$chr}{$start}};
			my $start_c = &convert($start);
			my $end_c = &convert($end);
			my $flag_no = 0;
			for(my $i = $start_c; $i <= $end_c; $i += 30){
				if(exists $mask{$chr}{$i}){
					if($mask{$chr}{$i} > 1){
						$flag_no = 1;
						last;
					}
				}
			}
			if($flag_no == 0){
				if($direction eq '+'){
					$BED_far{$chr}{$start} = $direction;
				}elsif($direction eq '-'){
					$BED_far{$chr}{$end} = $direction;
				}
			}
		}
	}
}



#---------------------------------------
# calculate enrichment around 3000bp
#---------------------------------------
my %PEAKS;
foreach my $chr(keys %BED_far){
	foreach my $pos(keys %{$BED_far{$chr}}){
		my $direction = $BED_far{$chr}{$pos};
		if($direction eq '+'){
			for(my $i = $pos - 3000; $i <= $pos + 3000; $i += 30){
				my $category = $i - $pos;
				push @{$PEAKS{$category}}, &GetWig($chr, $i);
			}
		}elsif($direction eq '-'){
			for(my $i = $pos + 3000; $i >= $pos - 3000; $i -= 30){
				my $category = $i - $pos;
				push @{$PEAKS{$category}}, &GetWig($chr, $i);
			}
		}
	}
}

#---------------------------------------
# output 上下5%ずつのデータは除去して平均を計算する
#---------------------------------------
print "position\taverage score\n";
foreach my $pos(sort {$a <=> $b} keys %PEAKS){
	my @sorted = sort {$a <=> $b} @{$PEAKS{$pos}};
	my $SUM = 0;
	my $NUM = 0;

	my $MAX_index = @sorted - 1;
	my $MIN_index = 0;
	# もしデータ数が10以上あるようなら、上位5%と下位5%を除く
	if(@sorted > 10){
		$MAX_index = int(@sorted * 0.95);
		$MIN_index = int(@sorted * 0.05);
	}

	for(my $i = $MIN_index; $i <= $MAX_index; $i++){
		$SUM += $sorted[$i];
		$NUM++;
	}
	printf "%d\t%f\n", $pos, $SUM / $NUM;
}



sub GetWig{
	my ($chr, $pos) = @_;
	for(my $i = 0; $i < 100; $i++){
		my $loc1 = $pos + $i;
		if(exists $WIG{$chr}{$loc1}){
			return $WIG{$chr}{$loc1};
		}
		my $loc2 = $pos - $i;
		if(exists $WIG{$chr}{$loc2}){
			return $WIG{$chr}{$loc2};
		}
	}
	return 0;
}



sub convert{
	my ($num) = @_;
	my $conv = int($num / 30) * 30 + 20;
	return $conv;
}


